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' Abstract. An embedding scheme is developed for the Dirac Hamiltonian H. 

Dividing space into regions I and II separated by surface S, an expression is derived for 
the expectation value of H which makes explicit reference to a trial function defined 
OO I in I alone, with all details of region II replaced by an effective potential acting on S 

and which is related to the Green function of region II. Stationary solutions provide 
^ , approximations to the eigenstates of H within I. The Green function for the embedded 

' Hamiltonian is equal to the Green function for the entire system in region I. Application 

of the method is illustrated for the problem of a hydrogen atom in a spherical cavity 
and an Au(001)/Ag/Au(001) sandwich structure using basis sets that satisfy kinetic 



o 



' balance. 
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1. Introduction 



X 

' There are many problems concerning electronic structure where attention is focussed on 
a small region of a larger system, at surfaces or defects in crystals being perhaps the most 
common. Let us call this region I, figure^ and the rest of the system region II. Although 
not of primary interest region II cannot be ignored, since in general the electron wave 
functions in I will be sensitive to the contents of region II. Some time ago Inglesfield 
|2j derived an embedding scheme which enables the single-particle Schrodinger equation 
to be solved explicitly only in region I. The influence of region II is taken into account 
exactly by adding an energy- dependent non-local potential to the Hamiltonian for region 
I, which constrains the solutions in I to match onto solutions in II. This embedding 
method has been developed into a powerful tool most notably for surface electronic 
structure problems [21 El IH IS] where it has found widespread application especially to 
situations where an accurate description of the spectrum of electron states is necessary. 
Examples include studies of image states [H], surface states at metals surfaces [Zj, static 
and dynamic screening [Hj, atomic adsorption and scattering at surfaces P, studies of 
surface optical response jTO] and field emission ^T]. Recent applications to transport 
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Figure 1. Schematic illustration of the embedded region I, the external region II, and 
the dividing surface S. 

problems have also been described fH]- For a review of the embedding method see 
Inglesfield [12]. 

In the case of materials containing heavier elements, relativistic effects can be 
significant jHj and lead to important deviations from the electronic structure as 
predicted by the Schrodinger equation - shifts in inner core levels of 5d elements 
are typically several 100 or 1000 eV, valence bands shifts are on the eV scale and 
spin orbit splitting is often measured in tenths of eV. Even ignoring the concomitant 
changes in electron wave functions these shifts can reorder levels and so affect calculated 
densities, fundamental to the determination of ground state properties within the density 
functional framework J3] . For this reason most of the conventional electronic structure 
techniques developed for accurately solving the single-particle Schrodinger equation in 
solids have subsequently been modified to deal with the Dirac equation, including the 
relativistic augmented plane wave method |161 [T7] , relativistic linear muffin-tin orbital 
method fH], relativistic augmented spherical wave method ^Hj and the relativistic 
multiple-scattering method [20], and each has subsequently been used in studying a 
diverse range of problems. The last method alone has formed the basis of calculations 
of photoemission [21], magnetocrystalline anisotropy [22]) hyperfine interactions [2S] and 
magnetotransport [21] amongst other topics. 

Inglesfield's embedding method has particular advantages that encourage its 
extension to the relativistic case. It permits the inclusion of extended substrates for 
surface and interface calculations, enables the study of isolated point defects in solids 
and being a basis set technique is highly flexible and permits full-potential studies with 
relative ease. At surfaces extended substrates (as against the use of the supercell or 
thin-film approximation in which the crystal is approximated by a small number of 
layers, typically 5-7) enable the proper distinction between surface states, resonances 
and the continuum of bulk states [2]. The behaviour of the W(llO) surface [2S] where 
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the addition of half a monolayer of Li is observed to increase the spin-orbit splitting of 
a surface state by ~ 0.5 eV (resulting in Fermi surface crossings separated by ~ 20% of 
the Brillouin zone dimension) typifies a type of problem a relativistic embedding scheme 
could address. Indeed each of the topics mentioned at the end of the previous paragraph 
are relevant at surfaces and/or interfaces, and could be usefully investigated within a 
relativistic embedding framework. 

In this paper we develop an embedding scheme for the Dirac equation that parallels 
Inglesfield's scheme for the Schrodinger equation. Inglesfield's starting point is the 
expectation value of the Hamiltonian using a trial wave function which is continuous 
in amplitude but discontinuous in derivative across the surface S separating I and II. 
The first order nature of the Dirac equation precludes the use of a similar trial function. 
Instead, in the following section we use a trial function in which the large component 
is continuous and the small component discontinuous across S. Continuity in the small 
component is restored when the resulting equations are solved exactly. Using the Green 
function for region II we are able to derive an expression for the expectation value purely 
in terms of the trial function in I. In section|21the application of the method is illustrated 
by calculating the eigenstates of a hydrogen atom within a cavity and in section |3] we 
determine the Green function for the embedded region. Section El briefly illustrates 
the method applied to a sandwich structure where relativistic effects are marked. We 
conclude with a brief summary and discussion. 

2. Embedding scheme 

In this section we consider region I joined onto region II (figure Q), and derive a 
variational principle for a trial wave function ip defined explicitly only within region 
I. We are primarily interested in the positive energy solutions of the Dirac equation |26j . 
and so we refer to the upper and lower spinors of the Dirac bi-spinor solutions as the 
large and small components of the wave function respectively. We notionally extend 
into II as Xy cin exact solution of the Dirac equation at some energy w, with the large 
components of ip and x (v^i Xi) niatching on the surface S separating I and II, but 
with no constraint upon the small components {ips and Xs), figure 121 The expectation 
value for the energy W is then 



where H = ccn-'p+ (3mc^ +V . (For clarity we omit the interaction oc (3(t B which appears 
in the relativistic density functional theory |^ neglecting orbital and displacement 
currents, where S is a "spin-only" effective magnetic field containing an external and 
exchange- correlation contribution. Its inclusion has no consequences for the derivation.) 
The first two terms in the numerator are the expectation value of the Hamiltonian 
through regions I and II, and the third the contribution due to the discontinuity in the 
small component of the wave function on S (in this and the following, surface normals 
are directed from I to II). 



W 



{^\H\^)i + {x\H\x)ii + ic^ js • [V^s - Xs] 
{"fWh + (xlx)ii 



(1) 
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Figure 2. The large component of the trial function is defined to be continuous across 
the surface S dividing regions I and II, but the small component can be discontinuous. 



We eliminate reference to x by introducing two relations. Firstly, for r G II, x 
satisfies the Dirac equation at energy w 

- ichct ■ Vx + [(3mc^ + V -w]x = (2) 
and differentiating with respect to w the energy derivative of x, X = dx/dw, satisfies 

- ichcy. ■ Vx + [Pmc^ + y — uj]x = x r G II. (3) 

Multiplying the Hermitian conjugate of the first equation by x from the right, 
multiplying the second from the left by subtracting and integrating over region 
II gives a relation between the normalisation of x iii II and the amplitude on S: 

{x\x)ii = ^ch drs-xWxs- (4) 
Js 

We have assumed that x vanishes sufficiently strongly at infinity. 

For the second relation we introduce the Green function (resolvant) G{r, r'; w) 
corresponding to equation 

- icha ■ VG + [pmc^ + V-w]G = -6{r - r') r, r' G II. (5) 

Multiplying the Hermitian conjugate of this equation by x from the right, and 
subtracting times equation Q, integrating over region II and then using the 
reciprocity of the Green function gives 

X{r) = ich / dr's ■ G{r, r^; w)a x(r^) r G II. (6) 
Js 
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We see that the Green function relates the amphtude of the wave function on S to the 
amphtude at any point within II. In particular, we can obtain a relation between the 
large and small components of x on S. Writing the 4x4 Green function as 

y G,i{r,r';w) Gss{r,r']w) J 

where each entry is a 2 x 2 matrix, substituting into equation (jH]), and rearranging the 
two equations coupling the small and large components of x gives 



where 



(rs) = ich ! dr's ■ T{rs, r'^; w)a xi(r^) (8) 
Js 



T{rs, r's; w) = G^rs, r's; w) + ich [ dr^ ■ G,,{rs, r'^; w)a- r(r^, r^; w) (9) 

Js 

It follows from (jH)) that the Green functions in (jH)) are the limiting forms of G{r, r^; w) 
as r ^ rs from within II. 

Equations (jH) and (jHl) are the desired results that enable us to express the 
expectation value IV in in terms of if alone. After substitution and use of the 
continuity of the large components xi = f\ on S we obtain 

{ip\H\ip)i + ich Jg drs ■ ^\cr ips — ich dr'g " |r — wt^ (jip\ 



W 



■ (10) 



{ip\ip)i - c^h^ drs ■ ip\(T dr'g ■ Ta ipi 
This is an expression for the expectation value of the energy W, given purely in terms 
of the trial function in region I and on the surface S, with all details of region II 
entering via F and its energy derivative. Following the convention in the non-relativistic 
embedding scheme we shall refer to F as the embedding potential. 

To see what this variational principle means in practice, we consider variations in 
(p\ whereby 

{S(f\H - W\ip)i + ich drs ■ Sipla - ich dr'g • |f + - w)f | a(fi 



6W 



{ip\ip)i - c^h^ drs ■ ^\cT dr's ■ Fcr (^i 



so that solutions if stationary with respect to arbitrary variations 5ip satisfy 

=Wi^ r G I (12a) 

ipsirs) = ichj dr's- ^^T{rs,r's;w) + iW -w)t{rs,r's;w)^aifii{r's). (125) 

The first expression indicates <y9 is a solution of the Dirac equation at energy W in region 
I. Comparing the second with (jHl) shows that ip also possesses the correct relationship 
between large and small components on S, the surface separating I and II, to match 
onto solutions in II. The term {W — w)r{w) provides a first order correction to T{w) so 
that the boundary condition is appropriate for energy W. 

In practice expression (fTUI) may be used to obtain solutions of the Dirac Hamiltonian 
by inserting a suitably parameterised trial function and varying the parameters to obtain 
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a stationary solution. This is conveniently achieved by expanding the trial solution in 
a finite basis of separate large and small component spinors 



Ni 



n=l 







iVs 



n=l 








' V'i(r-) 




ai 






VsW _ 




as 



(13) 



The matrix in the final expression is 4 by A^i + A^s, and the column vector contains 
the A^i + Ns coefficients. Substituting into (|Tn|l we find states that are stationary 
with respect to variations in the expansion coefficients {a\^n, ^s.n} are then given by the 
eigenstates of a generalised eigenvalue problem of the form 



where 



[Oll]nn' ^ 



H\\ His 




a\ 


= W 


On 




a\ 




Hsi Hss 




as 


Oss 




as 




nir) {Vir) + 






' fdrs- 4j} 




[ dr's- 


T{rs,r's;w) 


— wV 


{rs,r's]w) 


Js 




Js 







(14) 



' i)\^{r)ca- ■ p'ips,n'ir)dr + ich I dvg ■ ipl^{rs)(Tips,n'irs) 
I ' Js ' 

'^ln{r)c(T ■pilJi,n'{r)dr 
i^lnif) {Vir) - mc^) i)s,n'{,r)dr 
i^ini'^)i's,n'ir)dr 



-"^h^ j drs ■ V^/„(r5)cr / dr^ ■ r(rs, r^; w)crz/>i^„/(r^ 

^I,n(^)^s,n'(^)dT'- 



[Os 



(156) 
(15c) 
(15d) 

(15e) 
(15/) 



Of course the spectrum of the Dirac Hamiltonian is unbounded below, and care must 
be taken to prevent solutions collapsing to negative energies. This can be avoided 
through the use of a kinetically balanced basis j2Zj in which there is a one-to-one 
relationship between large and small component spinors, Ns = N\ = N, and where 
the small component spinors are given by 

^ps,ni'^) = cr ■ p ipi^r). (16) 

The upper half of the spectrum of the 2N eigenstates of ()14|1 then provide 
approximations to the spectrum of electronic states. 



3. Model application 



To illustrate the application of the relativistic embedding scheme we consider a model 
problem of a hydrogen atom within a spherical cavity, finding bound states of the Dirac 
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V(r) 




^ oo 



Figure 3. The model potential used to illustrate the relativistic embedding scheme. 



equation corresponding to the potential illustrated in figure El 
Vir) = 



-\/r 



r <R 
r > R 



where A = e^/(47reo) and Vq > 0. We choose this model as the bound states may also 
be found straightforwardly by alternative methods. Region I, the region to be treated 
explicitly, is the sphere of radius R centered on r = 0. The external region II where 
V{f) = Vq is replaced by an embedding potential acting on the surface of the sphere. 
The value of the embedding potential is most readily evaluated from equation (jSJ. A 
general solution to the Dirac equation at some energy w in region II and satisfying the 
appropriate boundary conditions is [201 



X[r) 



Xi[r) 
Xs{r) 



OA 



TT 



2kr \ -i-f^Ki^i{kr)Qj^{r) 



(17) 



where A = (k, /u), A = (— /t, /i), Qa^v) a spin-angular function, K^_^_i{z) a modified 
spherical Bessel function of the third kind f2Hl, chk = ^Jvnu^c^ — {w — Vq)^, 7^ = 
hck/ {w — Vq + mc^) and 



K > 

K < 



K 



:i8) 



The spherical symmetry of region II means the the embedding potential F may be 
expanded on 5* as 

r(r5, r's; w) = J2 rM^Airs)^\ir's) (19) 
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and substituting (fT^ and (fTTj) into (jSj) leads to 

"^^^ chR' K^^._{kR)' 

Using with the Green function for constant potential (w < Vo,r > r') 
{w — Vq + mc^) 



(20) 



G(r, r'; w) 



X 



/,+ i(fcr')r^A(r') 
i'lJl^i{kr')nx{r') 



(21) 



i^,+ l(A:r)^]A(r) 
-z7«,ir^-+i(A;r)fiA(^) 

gives the same result but after rather more involved manipulations. 1 {z) is a modified 
spherical Bessel function of the first kind. 

Because of the spherical symmetry we can determine separately states with a given 
angular character A. Using as a basis set for the large component spinors 



so that the small component spinors ensuring kinetic balance are 

Cn(^) = -fn.{r)n^{r), fn.{r) = h[{n + K)r^-' - r^] e 



r 

the matrix elements become 

R 



(22) 



(23) 



(A) 
11 







9n[r 



A 2 

— h mc 
r 



+h'c'R'gn{R)gn'{R) 



] nn' 



R 



fnK{r) 



X 



mc 



gn'ir)dr 
fn'K{r)dr 



R 



he / gn{r) 



(A) 



^11 



R 



dr 

dgrAr) 
dr 



r 



-fn'K{r) 



dr + hcgn{R)fn'K{R) 



+ -9n'[r) 
r 



dr 



gn{.r)gn'{r)dr -h c R gn{R)gn'{.R)^ k{w) 



nn' 



R 



fnK{r)fn'K{r)dr 



(24a) 
(246) 

(24c) 

(24 d) 

(24e) 



The eigenvalues only depend upon the quantum number k. In table ^ the lowest two 
eigenvalues of k = —1 symmetry (corresponding to the lsi/2 and 2si/2 of free hydrogen) 
are shown as a function of basis set size and for different values of the energy w at which 
the embedding potential is evaluated, for the case i? = 3, Vq = 10. For comparison also 
given are the values found by matching the external solution (|17|) to the regular internal 
solution, which can be expressed in terms of confluent hypergeometric functions |26j. For 
a given fixed w the eigenvalues converge from above to values that are equal or above the 
exact values. The further w lies from the eigenvalue, the larger the difference between 
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Table 1. Lowest two electron-like eigenvalues E — W — mc^ of symmetry of a 
hydrogen atom confined to a spherical cavity with radius R = 3 and confining potential 
V = 10, obtained using embedding potentials at different trial energies w. w = W 
indicates the solution has been iterated to ensure w coincides with the eigenvalue. The 
exact eigenvalues are those found by matching internal and external solutions at R. 
We use atomic units — h = m = 1 and c = 137.035 999 76 so that the corresponding 
free atom eigenvalues are -0.500 0067, -0.125 0021 





w = mf? — 0.5 


w = m(? 




2 
4 
6 
8 

'exact" 


-0.4111620, 1.698 0995 
-0.4451482, 0.912 9418 
-0.445 5519, 0.8914789 
-0.445 5532, 0.8912708 
-0.445 5532, 0.8908194 


-0.4111527, 1.6949300 
-0.4451439, 0.912 6817 
-0.445 5477, 0.8912219 
-0.445 5488, 0.8910141 
-0.445 5532, 0.890 8194 


-0.4111624, 1.689 6482 
-0.4396204, 0.9714775 
-0.445 5520, 0.8910268 
-0.445 5532, 0.890 8194 
-0.445 5532, 0.890 8194 



the limiting value for large basis sets and the correct value. However, the influence of 
the r terms in means the error is relatively small. When w = mc^, the lowest 
eigenvalue found with A'l = 8 is —0.445 5488 Ha and in error by only 0.000 0044 Ha, a 
factor 10^ smaller than the error in w. 

Differentiating (fTUj) with respect to the trial energy w shows the expectation value 
is stationary at w = W. In this case W is given by the solutions of 

_ {^\H\ip)i + ichjgdrs ■ [jfs - ichj^dr's ■ T{W)(Tipi] 

Eigenf unctions ip solving this equation satisfy the Dirac equation within I and the 
relationship between small and large components on S ()12&|) is exact. The final column 
in table HI shows the lowest two positive energy eigenvalues of (EKjl . again as a function 
of basis set size. The eigenvalues again converge from above, and by A'^i = 8 reproduce 
the exact values by at least 7 significant figures. It is worth noting that with this 
particular basis set increasing A^i much further leads to some numerical difficulties due 
to overcompleteness. For more accurate work a more suitable basis set should be used. It 
should also be noted that conventional finite basis set calculations using a basis satisfying 
kinetic balance can given eigenvalues that lie below exact limiting values by an amount 
of order jSZj, and similar behaviour is expected in this embedding scheme. 

4. Green function 

Most practical applications of the Schrodinger embedding scheme have actually used 
the Green function of the embedded system. This is a more convenient quantity when 
dealing with systems where the spectrum is continuous, such as at surfaces or defects 
in solids. We therefore consider the Green function for the embedded Dirac system. 
Differentiating (|T0|) with respect to w shows W is stationary when w = W, as would 
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be expected. In this case stationary solutions satisfy the embedded Dirac equation 

HyD{r) + j^dr'A{r,r';W)^{r') -W<^{r) = rel (26) 

where, introducing asirs), the component of cr in the direction normal to the surface 
5* (from I to II) at r^, the additional term A enforcing the embedding is 

A(r,r';w)= / drs5{r-rs) / dr'g5{r' - r'^) 
Js Js 



c^h'^asirs)T{rs, r'^; w)asir's) ichasirs)Sirs - r'g) 




(27) 



The corresponding Green function satisfies 

HG{r, r'- ^) + dr" A(r, r"; W)G{r", r'; W) -WG{r, r'; W) = -5{r - r') (28) 

for r, r' G /. A similar line of argument to that given by Inglesfield for the embedded 
Schrodinger equation shows that this Green function is identical for r,r' G / to the 
Green functions Gi+n for the entire system I+II. For simplicity assuming I+II constitute 
a finite system so that the spectrum is discrete, the Green function Gi+n is given by 

n " 

where W^^"*"^^ is the eigenvalue corresponding to eigenstate "^nir) of the entire system, 
normalised to unity over I+II. For a given W, the Green function solving (|28|) can 
be expanded in terms of the eigenstates (pnif] W) of the corresponding homogeneous 
equation 

Hifrrir; W) + dr'A(r, r'; W)^n{r'; W) - Wn{W)^n{r; W) = Q rel (30) 
normalised to unity over I, as 

(31) 

uyr,r,vv) ^ _ ^^^^^ . ^-^i; 

Clearly G has poles at IV = iy„(iy). At these energies (j30|) becomes the exact embedded 
Dirac equation (|26p so as we have seen the poles will occur at eigenstates of the entire 
system and the spectrum of G and Gi+n coincide. It remains to show the poles of G 
have the appropriate weight. The residue of G at Wn is 

The second term in the denominator is precisely the additional factor necessary to 
correctly normalise the states (see (jl)), (jHl)) so that 

l-idW^{W)/dW)\^^^ ^n(^)^n(r). [66) 

The residues of the Green function of the embedded system and those of the entire 
system are identical. Hence the two Green functions are identical for r, r' G I. 
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For practical calculations the Green function can be expanded using a double-basis 
of separate large and small component spinors: 



G{r,r';W) 











G{W) 









Air') 



t 



(34) 



The matrix elements of the matrix of coefficients G may be found by substituting into 
p8|l . multiplying from the right by the vector of basis functions, multiplying from the 
left by the Hermitian transpose of the vector of basis functions, and integrating over 
region I. This leads to 

WOn - Hn -H^, 



G{W) 



wo. 



(35) 



where the overlap and Hamiltonian matrices have their previous definitions 
with r = 0. 

As an illustration we calculate the local density of states for the confined hydrogen 
model at energies above Vq where the spectrum is continuous. Integrating over the 
embedded region this is given by 



n{W) = Im Tr G{W + iO^)0. 



IT 



(36) 



Figure m shows the Si/2 wave local density of states for R = 3, V = 1, calculated with 
varying number of basis functions. The basis functions (j22I), are not particularly 
appropriate for representing the continuum wave solutions, and so convergence is only 
achieved using a relatively large set; however the results serve to illustrate the systematic 
improvement that accompanies an increasing number of basis functions. The local 
density of states shows two resonances, the precursors of bound states that exist when 
any of R, Z or Vq are increased sufficiently. 



5. Application to an embedded monolayer 

As a further example, one that provides a test of the relativistic embedding scheme 
when applied to a more challenging problem, we use it to calculate the local density of 
states on a silver monolayer in a Au(001)/Ag/Au(001) sandwich structure. Using the 
embedding scheme only the region occupied by the Ag monolayer is explicitly treated. 
This is region I, with the two Au halfspaces to either side entering the calculation via 
embedding potentials expanded on planar surfaces. Then, using Bloch's theorem the 
calculation is performed within a unit cell containing one atom. The full technical 
details will be described elsewhere, but briefiy the Green function at two-dimensional 
wave vector K is expanded in a set of linearised augmented relativistic plane waves. 
We use large component basis functions 

{ipcr exp(i(l^ + G) ■ r) r G interstitial 

J2 [^^"u,{r) + B^'^ii^ir)] n^{r) r G muffin-tin (37) 



An embedding scheme for the Dirac equation 



12 




Figure 4. The local density of states as a function of energy E = W — mc^, 

integrated through a sphere of radius i? = 3 for the model system of a confined 
hydrogen atom with model parameters R = 3, Vq = 1. Different curves have been 
calculated with basis sets corresponding to the indicated number N\ of basis functions. 



where (pa- is a Pauli spinor, G = g + GzZ, with g a two-dimensional reciprocal lattice 
vector and Gz = nx2TT/D,n = 0, ±1, ±2, . . . and where D exceeds the width of the 
embedded region ensuring variational freedom in the basis. The function is the large 
component of the wavefunction that satisfies the radial Dirac equation for the spherically 
symmetric component of the potential at some pivot energy; is the energy derivative 
of Mk. The matching coefficients A^, ensure continuity of the basis function in 
amplitude and derivative at the muffin-tin radius. The small component basis functions 
are chosen to satisfy kinetic balance. 

Overlap and Hamiltonian matrix elements follow directly from these basis functions. 
The embedding potential is obtained from (jHl) using the general expression for a 
wavefunction outside a surface at wave vector K. This gives for the embedding potential 
describing the left Au half space 

r(rs,r's 



W + rac^ 



J2 (5- + 5+i?+-)(n-i?+-)"V, 



gag' a' 



X exp(i(K + g) ■ rs) exp{-i{K + g') ■ r'g) ip^ ® Lp„> 



{31 
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Figure 5. Local density of states on a Ag monolayer in a Au(001)/Ag/Au(001) 

sandwich structure calculated using a relativistic scattering method ( ) and 

relativistic embedding (a). The insets shows the calculation geometry (left) and the 
local density of states obtained from the Schrodinger equation (right). 

with 

^tgv = d{^- K^) ^.'^9' (39) 

The reflection matrix is found using standard layer-scattering methods [221 ■ ^ 
similar approach may be used to obtain an embedding potential for the right half space, 
which unlike the non-relativistic case differs from that for the left half space. 

Figure El compares the local density of states calculated using the relativistic 
embedding technique for an embedded Ag monolayer using embedding potentials 
corresponding to Au(OOl) with that found for an Au(001)/Ag/Au(001) sandwich 
geometry using relativistic scattering theory [SHI- The same Au and Ag potentials has 
been used in each case, and the local density of states found within the same muffin-tin 
volume. Therefore the results obtained with the two methods should be comparable, 
and we find that they are indistinguishable. This confirms that the embedding potential 
p8|l imposes the correct variational constraint upon wave functions for the embedded 
Ag monolayer so that they replicate the behaviour of an extended Au(001)/Ag/Au(001) 
sandwich structure. The inset in figure El shows the local density of states in the non- 
relativistic limit (c — > oo\ indicating the significant relativistic effects on the electronic 
structure which are correctly reproduced with this Dirac-embedding scheme. 

6. Summary and Discussion 

We have outline above an embedding scheme for the Dirac equation. It enables the Dirac 
equation to be solved within a limited region I when this region forms part of a larger 
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system, I+II. Region II is replaced by an additional term added to the Hamiltonian for 
region I, and which acts on the surface S separating I and 11. The embedding scheme 
is derived using a trial function in which continuity in the small component across 5* 
is imposed variationally. Expanding the wave function in a basis set of separate large 
and small component spinors, the problem of variational collapse is avoided by using a 
basis satisfying kinetic balance. Calculating the spectrum of a confined hydrogen atom, 
the method is shown to be stable and converge to the exact eigenvalues. We have also 
derived the Green function for the embedded Hamiltonian and illustrated its use in the 
continuum regime of the same confined hydrogen system and an Au/Ag/Au sandwich 
structure. These are demonstration calculations - future applications are likely to be to 
defects and surfaces of materials containing heavier (typically 5d) elements, within the 
framework of density functional theory. 

It is worthwhile to discuss further the use of a trial function that is discontinuous in 
the small component, since such a wave function gives rise to a discontinuous probability 
density and so would normally be dismissed in quantum theory. In non-relativistic 
quantum mechanics discontinuous trial functions are not permitted, since they possess 
infinite energy. However the Dirac equation is first order in p and as we have seen a 
perfectly regular expectation value of H results. Exploiting this freedom, the embedding 
scheme outlined above leads to solutions that are continuous in both large and small 
component only when the embedding potential r[rs,r'g;w) is evaluated at the same 
energy w as the energy W that appears in the Dirac equation itself, for then the 
relationship between small and large components on S inside (equation 112 &|1 and outside 
(equation |HJ coincide, the large components matching by construction. This may be 
achieved for example via the iterative scheme used in connection with equation (j25|) 
and the final column of tabled or explicitly when determining the Green function as in 
Section 0] These are the methods in which the non-relativistic embedding scheme has 
been most widely used. 

When w and W do not coincide, the solutions obtained via this embedding scheme 
will retain small components that are discontinuous across S. This may be unacceptable 
for certain apphcations, but the solutions continue to be valid approximations at least in 
as much as they provide estimates of the energies of the solutions of the Dirac equation, 
and so could suffice e.g. for interpreting spectroscopic measurements. This embedding 
scheme places no greater emphasis on a discontinuity in the small amplitude at S than on 
an incorrect (but continuous) amplitude elsewhere within the embedded region. It aims 
merely to optimise the energy of the state, and will retain a discontinuity in the small 
component if in doing so it can better (in terms of energy) approximate the solution 
inside the embedded region. In the non-relativistic embedding scheme the discontinuous 
derivative of the trial function implies a probability current (and electric current) that is 
discontinuous across the embedding surface. This is similarly unphysical, yet numerous 
applications such as those cited above have demonstrated the utility and accuracy of the 
method. Indeed, there have been many applications in which this scheme has been used 
to determine currents and or transport properties, such as in relation to surface optical 
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response [TUj and electron transport in electron waveguides or through domain walls 
[T^ . The reason for the success of these calculations is that they employed schemes in 
which the embedding potential was evaluated at the correct energy, ensuring that the 
derivative of the wavefunction was continuous across S. In practise there have been few 
calculations using the non-relativistic embedding scheme in which the energies did not 
coincide. 

There are a number of aspects of the method which are worthy of further 
consideration. We started with a trial function in which by construction the large 
component was continuous and the small component discontinuous across the surface 
S dividing I and II. We could have reversed these conditions, leading to a similar 
embedded Dirac equation but with a modified embedding term. The particular choice 
was motivated by the wish to have a theory which behaves reasonably in the limit 
c — > oo when the small component becomes negligible - a discontinuous amplitude is 
not permissible in trial solutions to the Schrodinger equation. However, the behaviour 
of the alternative formulation should be investigated. Perhaps in connection with this 
there is the question of the spectrum of negative energy solutions, to which we have 
paid scant attention. 

Exploring the c — > oo limit it might be possible to identify how to embed a 
relativistic region I within a region II treated non-relativistically - a 5d overlayer 
on a simple metal substrate might be a physical system where such a treatment is 
appropriate. There could be benefits in terms of computational resources expended if 
the embedding potential could be determined within the framework of a non-relativistic 
calculation, and there might also be useful insights in terms of simple models. Finally, 
in terms of implementation for realistic systems, some of the novel schemes for deriving 
embedding potentials 13] could certainly be adapted to the relativistic case. It would 
also be worthwhile to consider whether it is possible to use a restricted electron-like 
basis, in which the large and small component spinors are combined. This is common 
practice in most relativistic electronic structure calculations for solids when using basis 
set techniques (e.g. ^j), and would result in significant computational efficiencies. 
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